Naval  Research  Laboratory 


Washington.  DC  20375-5000 


NRL  Memorandum  Report  6704 


AD-A226  552 


T\  'T*  if  /■> 

.  '  \  i  v> 

(?*:*  Vt-.E'ZTF 

if  •'M  "  >"’Vv. 

v4  SEP  1?  1990 f-J  $ 

r& 

sJ 


A  Study  of  Confined  Diffusion  Flames 


J.  L.  Bllzey,*  K.  J.  Laskey,**  and  E.  S.  Oran 

Laboratory  for  Computational  Physics  and  Fluid  Dynamics 

'Berkeley  Research  Associates 
Springfield,  VA 

"Grumman  Space  Station  Program  Support  Division 
Reston,  VA 


September  4,  1990 


Approved  for  public  release;  distribution  ul 


**  -Uli 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  Ho  0704-0188 


bubJtc  fvoorting  burden  'or  this  cotleaion  of  infortndtion  is  dstimatM  to  ,ver*q*  '  hour  per  resoonsd.  including  the  time  to r  reviewing  instructions,  searching  tinting  data  sources 
gathering  and  maintaining  the  data  needed,  and  comoteetng  and  reviewing  the  collection  of  information.  Sena  comments  regarding  this  burden  estimate  or  any  other  r street  of  this 
c oltect ion  of  information,  including  suggestions  for  reducing  this  burden  to  Aashington  Headquarters  Services.  Directorate  tor  information  Operations  and  Reports,  tits  Jefferson 
Davis  Highway.  Suite  1204.  Arlington.  V*  22202-4302.  and  to  the  Office  of  Management  and  Budget.  Paoerwera  Reduction  Protect  (0204-0 1M>.  Washington.  DC  20503 


1.  AGENCY  USE  ONLY  (Lee ve  blink) 


4.  TITLE  AND  SUBTITLE 


2.  REPORT  DATE 

1990  September  4 


3.  REPORT  TYPE  AND  OATES  COVERED 


S.  FUNDING  NUMBERS 


A  Study  of  Confined  Diffusion  Flames 


6.  AUTHOR(S) 


J.  L.  Ellzey,  K.  J.  Laskey  and  E.  S.  Oran 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADORESS(ES) 

Naval  Research  Laboratory 
Washington,  DC  20375-5000 


PE  -  61153N 
PR  -  0N280-071 
TA  -  0110943 
WU  -  44153000 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

NRL  Memorandum 
Report  6704 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  AODRESS(ES) 

Office  of  Naval  Research 
800  N.  Quincy  Street 
Arlington,  VA  22217-5000 


10.  SPONSORING  /MONITORING 
AGENCY  REPORT  NUMBER 


11.  SUPPLEMENTARY  NOTES 


*  Berkeley  Research  Associates,  Springfield,  VA 

**  Grumman  Space  Station  Program  Support  Division,  Reston,  VA 


12*.  DISTRIBUTION  /AVAILABILITY  STATE.  IT 


12b.  DISTRIBUTION  CODE 


Approved  for  public  release;  distribution 
unlimited. 


13.  ABSTRACT  (Maximum  200  words) 

A  numerical  simulation  of ,  an  axisymmetric  confined  diffusion  flame  formed  between  a 
H2  -  N2  jet  and  coflowing  air  at  30  cm/s  is^  presented  in  this  paper.  For  the  initial  computations, 
the  restrictions  of  the  Burke)Schumann  theory  are  imposed  and  the  results  of  the  computation  are 
compared  with  the  analytical  solution  for  flame  location.  For  both  the  underventilated  and  overven¬ 
tilated  flames,  the  results  of  the  computations  are  in  excellent  agreement  with  the  analytical  solution. 
However,  the  flame  behavior  becomes  more  complex  as  the  restructions  are  relaxed.  When  variable 
diffusion  coefficients  and  densities  are  included  in  the  calculation,  small  radial  velocities  are  induced 
and  the  flame  interface  is  slightly  distorted.  When  heat  release  is  included,  the  flame  is  shorter  and 
an  unsteady  mixing  region  forms  at  the  fuel-oxidizer  interface.  The  instabilities  are  damped  when 
viscous  effects  are  included.  Large-scale  instabilities  form  in  the  oxidizer  region  with  a  frequency  of 
approximately  15-20  Hz  when  gravity  is  included  in  the  calculation. 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

UNCLASSIFIED 


NSN  7540-01  -280-5500 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

UNCLASSIFIED 


19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

UNCLASSIFIED 


IS.  NUMBER  OF.  PAGES 

31 


16.  PRICE  CODE 


20.  LIMITATION  OF  ABSTRACT 


Standard  Form  298  (Rev  2-89) 

*»f«crio*<i  Dy  ANSI  Std 
298' 02 


CONTENTS 


I.  Introduction  .  1 

II.  Numerical  Methods  and  the  Model  Structure  .  2 

ED.  Application  of  the  Algorithm  to  a  Confined  Diffusion  Flame  .  11 

IV.  The  Burke-Schumann  Flame  .  12 

V.  Elimination  of  the  Burke-Schumann  Restrictions  .  14 

VI.  Conclusions .  16 

Acknowledgements  . 19 

References  .  20 


A.:-  -  7 

r-ns 

G ;  u 
u-  !.'■ 

Jt. 

^  i  i 

-f'3  U  I 

it.  ••  __  _  __i 

By 

I 

Dot-ib 

-  tio '  /  | 

A 

.'  ty  Codes 

t 

Av, . .  1  '<  ,  or 

Oist 

Si  'f_-CI.ll  j 

fl-l 

)  1 

1  ! 

!  i 

! 

iii 


A  STUDY  OF  CONFINED  DIFFUSION  FLAMES 


I.  Introduction 

The  analytical  work  of  Burke  and  Schumann  (1)  has  formed  many  of 
our  fundamental  ideas  about  laminar  diffusion  flames.  In  the  original  Burke- 
Schumann  problem,  axisymmetric  coflowing  streams  of  fuel  and  oxidizer  flow 
through  a  confined  duct,  and  the  velocities,  densities,  and  diffusion  coeffi¬ 
cients  of  the  fuel  and  oxidizer  Eire  equal.  Reaction  is  instantaneous,  resulting 
in  a  flame  sheet  of  infinitesimal  thickness  in  which  the  reaction  rate  is  ef¬ 
fectively  controlled  by  the  diffusion  rate.  The  solution  to  this  problem  is  an 
equation  which  can  be  solved  for  the  location  of  the  flame  front.  In  the  orig¬ 
inal  analysis,  Burke  and  Schumann  chose  the  diffusion  coefficients  in  order 
to  obtain  good  agreement  with  experiments.  Later  work  (2)  extended  the 
theory  to  describe  the  flame  interface  for  unequal  velocities  and  diffusion  co¬ 
efficients.  Further  extensions  of  the  theory  predict  the  behavior  of  multiple, 
coupled  diffusion  flames  (3,4). 

Many  of  the  restrictive  assumptions  in  the  Burke- Schumann  analysis 
can  be  removed  by  a  numerical  solution  of  the  reactive  flow  equations.  In 
particular,  there  are  a  number  of  steady-state  numerical  solutions  that  sim¬ 
ulate  laminar  diffusion  flames.  Gosman  et  al.  (5)  solved  the  two-dimensional 
steady-state  equations  for  a  case  in  which  all  of  the  diffusion  coefficients  were 
the  same  and  the  Lewis  number  was  unity.  Mitchell  et  al.  (6)  numerically 
simulated  steady-state  flames  with  nonunity  Lewis  numbers  but  kept  the 
basic  idea  of  the  flame  sheet  model. 

This  paper  describes  a  time-dependent,  axisymmetric,  compressible  nu¬ 
merical  model  which  is  designed  specifically  to  simulate  the  nonsteady  be¬ 
havior  of  diffusion  flames.  It  contains  submodels  for  finite-rate  chemistry, 
viscosity,  thermal  conduction,  and  the  temperature  dependence  of  material 
properties  such  as  specific  heats  and  diffusion  coefficients.  As  one  of  the 
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first  uses  of  the  diffusion  flame  model,  we  simulate  &  Burke- Schumann  flame 
and  remove  the  restrictions  individually.  We  present  results  for  a  classic 
Burke- Schumann  flame  with  all  of  the  restrictions  included  in  the  analy¬ 
sis  and  compare  to  the  anlytical  solution.  Then  we  include  the  following 
sequentially: 

1.  Variable  density  and  diffusion  coefficients,  radial  convection,  and  axial 
diffusion, 

2.  Heat  release, 

3.  Viscous  effects, 

4.  Gravitational  effects. 

This  set  of  computations  is  a  benchmark  of  the  model  that  is  currently  being 
applied  to  more  complex  transitional  diffusion  flames.  Although  practical 
examples  of  Burke-Schumann  flames  are  not  as  abundant  els  turbulent  or 
unsteady  flames,  they  represent  an  important  class  of  problems  which  can 
be  studied  through  theory,  computation,  and  experiment. 

Besides  the  applications  of  the  model  to  the  Burke-Schumann  problem, 
this  paper  also  describes  the  numerical  model  in  detail.  Specific  notable  fea¬ 
tures  of  the  model  axe  its  time  dependence,  the  finite-rate  chemical  model,  the 
temperature  dependence  of  the  transport  coefficients,  and  the  nonunity  Lewis 
number.  The  new  elements  of  the  numerical  model  that  make  such  compu¬ 
tations  possible  are  the  BIC-FCT  algorithm  used  to  compute  the  convection 
and  the  parametric  diffusion-reaction  (PDR)  model  used  for  the  finite-rate 
chemistry.  These  features  and  algorithms  are  described  in  some  detail  in  the 
next  section  of  this  paper. 

II.  Numerical  Methods  and  the  Model  Structure 

The  numerical  model  used  for  this  study  is  based  on  those  originally 
developed  by  Patnaik  et  al.  (7)  to  simulate  low-speed  premixed  flames  and 
by  Laskey  (8)  to  simulate  jet  flames  .  The  program  solves  the  equations 
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de  =  pcv  dT  .  (6) 

The  various  quantites  used  in  these  equations  and  throughout  this  article  are 
defined  in  the  accompanying  Nomenclature. 

Equations  (1)  -  (4)  contain  terms  representing  convection,  thermal  con¬ 
duction,  species  diffusion,  chemical  reactions,  and  viscosity.  These  equations 
are  then  rewritten  in  terms  of  finite-difference  approximations  on  an  Eulerian 
mesh  and  solved  numerically  for  specified  boundary  and  initial  conditions. 
The  accuracy  of  the  solution  is  determined  by  the  specific  finite-difference 
algorithm,  the  spatial  resolution  set  by  the  computational  grid,  and  the  tem¬ 
poral  resolution  set  by  the  timestep.  There  is  a  wide  range  of  important 
spatial  and  temporal  scales  in  reacting  flow  problems.  Because  it  is  not  usu¬ 
ally  possible  to  resolve  phenomena  on  all  of  these  scales,  the  smallest  scales 
must  be  modeled  phenomenologically. 

However,  as  computational  capabilities  and  model  inputs  improve,  it 
should  be  possible  to  replace  certain  submodels  by  more  accurate  or  faster 
submodels.  For  example,  assuming  that  rate  of  diffusion  is  much  less  than 
the  reaction  rate  helps  justify  using  a  global  reaction  mechanism.  Using 
such  a  global  reaction  is,  however,  approximate  and  we  believe  that  within  a 
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few  years  it  will  be  possible  to  include  a  detailed  set  of  elementary  reaction 
rates.  Therefore,  the  computer  program  is  designed  in  a  modular  form  so  that 
particular  submodels  can  be  updated  in  a  relatively  straightforward  way.  In 
the  computer  program,  algorithms  representing  different  physical  processes 
are  solved  separately  and  then  the  results  are  combined,  as  summarized 
below.  More  detailed  descriptions  are  presented  by  Laskey  (8). 

Convection 

The  solution  to  the  convective  terms  in  Eqs.  (1)  -  (4)  is  obtained  us¬ 
ing  the  new  algorithm,  Barely  Implicit  Correction  to  Flux-Corrected  Trans¬ 
port  (BIC-FCT)  that  was  developed  to  solve  the  convection  equations  for 
low-velocity  flows  (9).  The  Flux- Corrected  Transport  algorithm  itself  is  an 
explicit,  finite-difference  algorithm  that  is  constructed  to  have  fourth-order 
phase  accuracy  (10).  Through  a  two-step  predictor-corrector  algorithm,  FCT 
ensures  that  all  conserved  quantities  remain  monotone  and  positive.  The 
FCT  procedure  is  to  first  modify  the  properties  of  a  high-order  algorithm 
by  adding  diffusion  during  a  convection  step  and  then  to  subtract  out  the 
diffusion  in  an  antidiffusion  phase.  In  addition,  fluxes  are  limited  to  ensure 
that  no  new  unphysical  maxima  or  minima  are  added  during  the  convection 
process. 

However,  because  FCT  is  an  explicit  algorithm,  the  numerical  timestep 
required  for  accuracy  and  stability  is  limited  by  the  velocity  of  sound  accord¬ 
ing  to  the  Courant-Friedrichs-Lewy  condition,  At  <  min(Ax/cs).  To  avoid 
this  restriction,  which  would  make  computations  of  slowly  evolving  flows  pro¬ 
hibitively  expensive,  the  convection  equations  are  usually  solved  implicitly. 
This  filters  out  the  sound  waves  from  the  equation  and,  therefore,  removes 
the  sound-speed  condition.  Patnaik  et  al.  (7)  developed  BIC-FCT  so  that 
the  timestep  would  be  limited  by  the  fluid  velocity  and  not  the  sound  speed. 
This  implementation  has  great  advantages  for  computations  of  slow  flows 
because  one  BIC-FCT  timestep  costs  the  same  as  one  regular  FCT  explicit 
timestep,  but  the  size  of  the  timestep  might  be  a  factor  of  50  to  100  times 
greater. 
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BIC-FCT  is  based  on  the  idea  proposed  by  Casulii  and  Greenspan  (11) 
that  only  the  terms  containing  the  pressure  in  the  momentum  equation  and 
the  velocity  in  the  energy  equation  must  be  treated  implicitly  in  order  to 
avoid  the  sound-speed  limitation  on  the  timestep.  BIC-FCT  has  three  steps. 
In  the  first  step,  the  conservation  equations  are  solved  explicitly  with  FCT 
using  a  relatively  large  timestep  governed  by  fluid  velocity.  In  the  second 
step,  the  energy  and  momentum  equations  are  rewritten  in  terms  of  a  pres¬ 
sure  correction,  Sp.  These  equations  can  be  manipulated  such  that  only  one 

elliptic  equation  for  Sp  must  be  solved.  In  the  third  step,  final  values  of 
momenta  and  energy  are  obtained  by  adding  the  pressure  correction  terms. 

Pafnaik  et  al.  (7)  have  used  this  algorithm  in  a  two-dimensional  flame  pro¬ 
gram  to  investigate  laminar  instabilities  in  premixed  flames. 

The  form  used  for  the  relation  between  the  change  in  internal  energy 
and  the  change  in  the  pressure,  required  in  the  second  step  of  BIC-FCT, 
can  be  derived  from  Eq.  6  and  simplifies  under  certain  assumptions  (8).  For 
example,  if  all  of  the  constituents  have  the  same  temperature  dependence, 
we  can  write 

*-  (7) 

where  7 m  is  a  mixture  quantity.  The  second  term  in  Equation  (7)  is  depen¬ 
dent  on  the  change  in  the  local  species  concentration  and  is  important  only 
in  the  reaction  zone.  In  our  simplified  reaction  submodel,  we  do  not  include 
the  intermediate  species.  Consequently,  this  term  cannot  be  represented  ac¬ 
curately  and  we  do  not  include  it  explicitly.  In  the  reaction  zone,  we  are 
primarily  interested  in  representing  a  finite  flame  thickness  and  reproducing 
the  correct  maximum  temperature.  We  calibrate  the  reaction  submodel  to 
predict  a  realistic  flame  temperature  and  flame  thickness  without  this  term 
included. 

Molecular  Diffusion 

An  algorithm  for  molecular  diffusion  has  been  formulated  to  estimate 
the  molecular  diffusion  fluxes  without  having  to  solve  a  full  matrix  problem. 
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The  change  in  species  concentration  for  each  species  k  due  to  molecular 
diffusion  is 


dnk 

dt 


-V  •  nkvk 


(8) 


where  the  diffusion  velocity,  v*,  is  calculated  from  Fick’s  Law, 


Vfc  =  ~^rDkmVXk  (9) 

A* 

and  then  corrected  by  a  procedure  described  by  Kee  et  al.  (12)  to  satisfy 
the  requirement  that  the  sum  of  the  diffusion  fluxes  is  zero.  This  method  is 
algebraically  equivalent  to  the  first  iteration  of  the  DFLUX  algorithm  (13), 
an  iterative  approach  that  solves  for  diffusion  velocities  to  optional  accuracy. 
The  change  in  total  energy  density  due  to  molecular  diffusion  alone  is 

dE  n‘p 

-q£  =  -V  •  ^2,  nkhkvk  ■  (10) 

*=i 


This  energy  term  is  calculated  during  the  diffusion  algorithm  but  is  added 
to  the  total  energy  at  the  end  of  the  time  step. 

The  explicit  finite-differencing  procedure  applied  to  this  term  introduces 
a  numerical  stability  condition, 

+  2^0  <  2  ’  (U) 

where  Dkm  is  the  diffusion  coefficient  for  species  k  diffusing  into  a  mixture. 
To  maintain  stability,  this  condition  may  require  a  timestep  smaller  than  that 
required  by  the  convection,  which  adds  substantially  to  the  cost  of  the  calcu¬ 
lation.  To  avoid  this  problem,  the  diffusion  term  is  evaluated  several  times 
during  a  convection  timestep.  This  is  especially  important  if  the  elevated 
temperature  of  the  reacting  flow  results  in  higher  diffusion  coefficients. 

Binary  diffusion  coefficients  are  calculated  from  kinetic  theory  and  axe 
in  the  following  form 

Dkt  =  TBkl  ,  (12) 

n 
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where  Aki  and  Bui  depend  on  species  k  and  /.  Values  for  Aki  and  Bki  have 
been  tabulated  by  Kailasanath  et  al.  (14).  The  diffusion  coefficient  of  species 
Ar  in  a  mixture  of  ntp  species  is  calculated  according  to 


Dk 


i-n 


kk,tk 


(13) 


where  Yk  is  the  mass  fraction  of  species  k,  Xk  is  the  mole  fraction  of  species 
k,  and  Dk,i  is  the  diffusion  coefficient  of  species  k  diffusing  into  species  /  (15). 

Thermal  Conduction 

A  two-dimensional  model  has  also  been  formulated  to  simulate  thermal 
conduction.  Restricting  our  attention  only  to  the  Fourier  conduction  term, 
the  energy  equation  appears  as 


dt 


-V  ■  (kVT)  . 


(14) 


As  with  the  molecular  diffusion  algorithm,  the  use  of  explicit  finite  differenc¬ 
ing  introduces  a  stability  limit  for  the  thermal  conduction  calculation, 


*A t  (  1  1  \  1 

pcp  \Ax2  +  Ay2/  <  2 


(15) 


where  */pcp  is  the  thermal  diffusivity.  The  thermal  conduction  term  is  eval¬ 
uated  several  times  during  each  convection  time  step  in  the  same  manner  as 
the  molecular  diffusion  term.  Thermal  conductivities,  «*,  for  the  individual 
species  were  calculated  from  kinetic  theory  over  the  temperature  range  300 
K  to  3300  K,  and  these  values  were  fit  to  a  third-order  polynomial.  The 
mixture  thermal  conductivity  is  then  calculated  using  the  expression  from 


Kee  et  al.  (12) 


*=  9 


1 

Y2xkKk  +  — 


k=l 


£ 
*= i 


(16) 
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Viscous  Stress 


The  momentum  transport  associated  with  viscous  diffusion  only  is 


dpv 

~dt 


=  —V  •  r 


(17) 


where 


(V-  v)I-fi  [(Vv)  +  (Vv)T 


(18) 


and  £  is  the  second  coefficient  of  viscosity  and  is  assumed  to  equal  zero. 
Equation  17  is  rewritten  in  terms  f  an  explicit  finite  difference  approxima¬ 
tion  which  introduces  a  numerical  stability  condition, 


/iAt  /  1  IN 
p  \Ai2  Ay2/ 


(19) 


The  values  for  pk  were  calculated  from  kinetic  theory  over  the  tempera¬ 
ture  range  300  K  to  3000  K  and  fit  to  a  third-order  polynomial.  The  mixture 
viscosity  was  calculated  from 


XkPk 


E 


where  $kj  is  a  weighting  factor 


and  Mk  is  the  molecular  weight  of  species  k  (16). 


(20) 


(21) 


Model  for  Chemical  Reactions 

Ideally,  we  would  like  to  simulate  the  chemical  reaction  by  including 
a  detailed  set  of  elementary  reactions  to  describe  the  production  of  the  in¬ 
dividual  species  and  the  energy  release  in  the  flame.  However,  the  cost  of 
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computer  time  and  memory  mak^s  this  prohibitive  for  problems  in  which 
the  flows  are  complex.  We  have,  therefore,  used  the  parametric  diffusion- 
reaction  (PDR)  model  that  reflects  many  of  the  major  characteristics  of  a 
more  detailed  calculation. 

In  the  original  flame-sheet  model  proposed  by  Burke  and  Schumann 
(1),  the  fuel  and  oxidizer  react  completely  and  are  not  permitted  to  coexist. 
Thus,  the  flame  is  an  infinitesimal  interface  between  regions  of  fuel  and 
oxidizer.  In  the  PDR  model  developed  by  Laskey  (8),  a  single  global  reaction 
mechanism,  2i?2  +  02  2i?20,  is  used  but  the  fuel  and  oxidizer  do  not  react 

instantaneously.  Instead,  the  reaction  occurs  over  a  finite  time  interval. 

In  a  real  flame,  the  increase  in  temperature  in  the  reaction  zone  is  a 
result  of  the  change  in  internal  energy  of  the  local  mixture  which  includes 
reactants,  product,  and  intermediates.  In  our  simulation,  we  do  not  include 
the  intermediates  and,  as  a  result,  the  specific  heat  of  the  gas  is  not  repre¬ 
sented  accurately  in  the  flame  zone.  Even  if  we  obtain  a  realistic  overall  rate 
for  the  global  reaction,  the  temperature  in  the  flame  zone  is  too  high  unless 
we  compensate  for  the  inaccuracy  in  the  specific  heat. 

The  two  important  physical  characteristics  of  the  flame  zone  are  the 
maximum  temperature  and  thickness.  We  calibrated  the  reaction  rate  in  a 
one-dimensional  transient  diffusion  flame  such  that  a  thin  reaction  zone  de¬ 
veloped  and  adjusted  the  heat  of  formation  such  that  the  maximum  temper- 
ture  was  approximately  the  adiabatic  flame  temperature  for  a  stoichiometric 
mixture  of  fuel  and  oxidizer. 

Coupling 

A  complete  solution  to  the  governing  equations  requires  solving  the 
terms  for  individual  processes  as  well  as  accounting  for  the  interaction  among 
the  processes.  In  the  calculations  presented  below,  we  use  timestep  splitting, 
which  assumes  that  the  net  effect  of  all  the  processes  is  the  sum  of  the  so¬ 
lutions  to  individual  processes.  This  technique  is  valid  if  the  changes  in  the 
dependent  variables  during  a  timestep  axe  small.  Table  1  is  an  outline  that 
shows  the  order  of  the  computations  during  one  timestep  in  the  computer 
code. 
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Table  1.  Outline  of  Diffusion  Flame  Code 


Initialize  Variables 
*  Increment  time 

1.  Thermal  Conduction 

Integrate  from  t  to  t  +  A t: 

Calculate  Aci 

Do  not  update  any  variables 
(Subcycle  as  necessary) 

2.  Ordinary  Diffusion 

Integrate  from  t  to  t  +  At: 

Only  update  {n,(x)} 

Calculate  At2 
(Subcycle  as  necessary) 

3.  Viscosity 

Integrate  from  t  to  t  +  At: 

Only  update  pv 
Calculate  A«3 

4.  Chemical  Reactions 

Integrate  from  t  to  t  +  At: 

Only  update  {n,(x)} 

Calculate  Ae* 

5.  Convective  Transport 

Integrate  from  t  to  t  +  At: 
x  direction  transport 

Update  p,  pv,  E ,  n* 
y  direction  transport 

Update  p,  pv,  E,  ri{ 

Implicit  correction-  update  p,  e,  and  E 
Start  New  Timestep  (go  to  *  above) 


Due  to  different  requirements  of  accuracy  and  stability,  the  type  of  cou¬ 
pling  used  for  lower  velocity  implicit  calculations  is  different  from  that  needed 
for  higher  velocity  explicit  calculations.  General  information  on  these  ap¬ 
proaches  is  described  in  some  detail  in  Reference  (10),  Chapter  13.  In  these 
implicit  computations,  the  changes  in  pressure  or  internal  energy  resulting 
from  the  individual  processes  should  not  be  added  into  the  solution  as  soon 
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as  they  are  computed,  but  instead  should  be  accumulated  over  the  timestep. 
The  entire  change  in  internal  energy  is  then  included  in  the  fluid  convection 
step.  The  specific  coupling  technique  used  in  this  program  (17)  allows  larger 
changes  in  variables  per  timestep  while  maintaining  numerical  stability. 

III.  Application  of  the  Algorithm  to  a  Confined  Diffusion  Flame 

The  numerical  procedure  described  above  was  used  to  simulate  several 
confined  diffusion  flames.  The  geometry  used  in  the  calculations  (Figure  1) 
consists  of  a  an  inner  jet  of  radius  a  and  an  outer  annular  region  between  the 
jet  and  the  walls  at  radius  b.  Typically,  fuel  flows  through  the  inner  jet  and 
oxidizer  flows  through  the  outer  annular  region.  The  appropriate  boundary 
conditions  for  this  geometry  are 

1.  r  =  0  is  a  line  of  symmetry, 

2.  r  =  b  is  a  solid,  adiabatic,  free-slip  wall, 

3.  z  =  0  is  an  inflow  boundary  where  the  concentrations,  velocities,  and 
temperatures  of  the  fuel,  oxidizer,  and  inert  are  specified, 

4.  z  =  Z\  is  an  outflow  boundary  where  the  pressure  is  adjusted  to  equal  1 
atmosphere. 

For  the  B urke- S chumann  flame,  a  —  1,  b  =  2,  z\  =10  cm.  The  com¬ 
putational  domain  consists  of  a  32  x  88  grid.  The  grid  spacing  is  uniform 
in  the  radial  and  the  axial  directions.  Calculations  on  finer  grids,  such  as 
64  x  176,  resulted  in  smoother  flame  interfaces  but  did  not  change  the  result 
significantly.  The  computational  time  step  is  1  ms. 

For  the  simulations  of  a  confined  diffusion  flame  without  all  of  the  Burke- 
S chumann  restrictions,  a  =  0.5,  b  =  2.5,  z\  —  10  cm.  The  grid  consisted  of 
64  x  88  cells  with  fine  cells  concentrated  around  the  jet  exit.  In  the  radial 
direction,  the  grid  spacing  is  approximately  0.02  cm  from  the  centerline  to 
r  =  0.7  cm  and  then  expands  gradually  to  a  grid  spacing  of  0.12  at  r  =  2.5 
cm.  In  the  axial  direction,  the  grid  spacing  is  uniform  through  the  domain 
and  equals  0.11  cm.  A  typical  timestep  is  10  fis. 
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Computational  Time  Requirements 

A  two-dimensional  simulation,  which  includes  convection,  chemical  reac¬ 
tion,  molecular  diffusion,  viscous  diffusion,  and  conduction,  requires  approx¬ 
imately  25  ps  per  grid  point  per  timestep  on  a  Cray  Y-MP.  The  convection 
algorithm  requires  approximately  twice  the  cpu  time  of  either  the  molecular 
diffusion  or  the  viscous  diffusion  algorithms  and  four  times  that  of  the  con¬ 
duction  algorithm.  The  parametric  diffusion-reaction  flame  model  requires 
insignificant  cpu  time. 

IV.  The  Burke-Schumann  Flame 

Burke  and  Schumann  (1)  found  a  solution  for  a  set  of  equations  that  give 
the  location  of  the  flame  interface  for  a  laminar  diffusion  flame  under  a  certain 
set  of  limiting  conditions.  The  Burke-Schumann  analysis  of  the  laminar 
diffusion  flame  and  a  comparison  of  their  analysis  to  the  computatuional 
results  are  presented  in  this  section. 

In  order  to  solve  the  equations,  Burke  and  Schumann  invoked  a  number 
of  simplifying  assumptions: 

1.  The  velocities  of  the  fuel  and  oxidizer  are  equal  and  uniform  everywhere, 

2.  The  radial  velocity  is  zero, 

3.  The  densities  and  diffusion  coefficients  are  equal  for  all  components, 

4.  Radial  diffusion  is  much  greater  than  axial  diffusion, 

5.  Reaction  takes  place  at  an  infinitesimal  flame  sheet. 

With  these  assumptions,  the  conservations  equations  can  be  reduced  to  a 
single  species  equation  which  is  solved  with  the  following  boundary  condi¬ 
tions: 

1.  r  =  b  is  a  solid  wall, 

2.  r  =  0  is  a  line  of  symmetry, 

3.  At  z  =  0,  the  compositions  of  the  fuel  and  oxidizer  streams  are  specified. 
The  analytic  solution  to  the  equation  yields  the  location  of  the  flame  surface 
as  a  function  of  a/6  and  the  initial  concentrations  of  fuel  and  oxidizer. 
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These  assumptions  enforce  various  unrealistic  restrictions  on  the  flow 
field.  If  the  velocity  is  uniform  across  the  radius  of  the  tube,  then  the  no¬ 
slip  condition  cannot  be  imposed  at  the  wall  boundary  and,  as  a  result,  a 
parabolic  velocity  profile  typical  of  confined  flows  cannot  develop.  The  as¬ 
sumption  of  equal  densities  requires  that  one  mole  of  fuel  reacts  with  s  moles 
of  oxidizer  to  form  1  +  a  moles  of  product.  Finally,  in  a  real  flame,  the  vol¬ 
umetric  expansion  associated  with  heat  release  distorts  the  one-dimensional 
flow  resulting  in  a  radial  component  to  the  velocity  and  a  nonuniform  density 
field.  Consequently,  the  heat  release  and  expansion  are  not  considered  in  the 
Burke-Schumann  analysis. 

Figure  1  shows  two  general  cases  that  can  be  solved  with  the  Burke- 
Schumann  approach:  the  underventilated  flame,  which  has  insufficient  oxi¬ 
dizer  for  complete  burning,  and  the  overventilated  flame,  which  has  excess 
oxidizer.  In  the  overventilated  case,  the  fuel  is  completely  consumed  in  the 
reaction  and  the  flame  surface  is  closed  at  the  centerline  of  the  jet.  In  the 
underventilated  case,  the  flame  surface  bends  outward  and  is  attached  to  the 
outer  wall.  The  two  different  cases  may  be  obtained  by  changing  either  the 
ratio  a/6  or  the  composition  of  the  fuel  or  oxidizer  stream. 

While  the  details  of  most  flames  cannot  be  represented  realistically  by 
the  results  of  the  Burke-Schumann  analysis,  the  predictions  of  the  analysis 
are  surprisingly  good  for  steady,  laminar  flames.  In  addition,  it  provides  an 
analytical  result  against  which  to  test  the  numerical  model. 

Simulation  of  the  Burke-Schumann  Flame 

For  the  first  cases  considered  in  this  study,  the  ratio  of  the  radii,  a/6,  is 
0.5  (a  =  1  cm,  6  =  2  cm),  and  the  fuel  flows  in  the  inside  jet  and  oxidizer  flows 
in  the  outside  annular  region.  The  velocities  of  the  fuel  and  oxidizer  streams 
are  uniform  and  equal  to  10  cm/s.  The  densities  of  the  two  streams  are  equal 
but  the  composition  was  varied  by  diluting  the  fuel  or  oxidizer  stream  with 
an  inert  gas.  All  diffusion  coefficients  were  equal  to  an  equivalent  mixture  of 
#2  and  N2  diffusing  into  O2. 
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Figure  2  shows  a  sequence  of  fuel  and  oxidizer  concentration  contours 
for  the  numerical  simulation  of  the  Burke-Schumann  flame  as  it  evolves  from 
the  initial  condition  to  a  steady  state.  At  t  =  0,  pure  fuel  exists  at  r  <  1 
cm  and  an  oxidizer  mixture,  consisting  of  one  part  oxidizer  and  one  part 
inert,  exists  at  r  >  1  cm.  Since  there  is  an  overall  excess  of  oxidizer,  the 
fuel  and  oxidizer  interface  moves  inward  during  time  steps  0  to  1000  because 
the  fuel  is  completely  consumed  in  the  reaction  at  the  flame  surface.  In  all 
cases,  fuel  and  oxidizer  do  not  coexist  because  the  Burke-Schumann  flame 
sheet  model  assumes  that  the  reaction  goes  to  completion.  At  time  step  1000 
which  equals  1  second  of  physical  time,  the  concentration  field  has  reached 
a  steady  state. 

In  Figure  3a,  the  computed  contours  for  1%  of  the  inlet  fuel  concen¬ 
tration  and  1%  of  the  inlet  oxidizer  concentration  are  superimposed  on  the 
Burke-Schumann  solution  for  the  flame  front.  The  analytic  solution  is  the 
interface  between  the  fuel  and  oxidizer  regions  and  should  correspond  to  the 
zeroth  contour  for  either  the  fuel  or  the  oxidizer,  and  it  should  occur  between 
the  two  1%  contours.  Figure  3b  shows  a  similar  calculation  for  a  fuel  mix¬ 
ture  of  one  part  fuel  and  three  parts  by  volume  of  inert  reacting  with  pure 
oxidizer.  In  both  cases,  the  computed  flame  front  is  within  one  cell  of  the 
analytical  solution. 

A  similar  computation  was  conducted  for  an  underventilated  Burke- 
Schumann  flame.  In  this  case,  the  inlet  jet  was  pure  fuel  and  the  oxidizer 
consisted  of  one  part  of  oxidizer  to  four  parts  inert.  The  steady  state  fuel 
contours  are  shown  in  Figure  3c  with  the  analytic  solution  superimposed.  In 
this  case,  the  flame  bends  outward  and  attaches  to  the  outer  wall.  Again, 
the  flame  shape  and  height  are  within  the  accuracy  of  the  calculation. 

V.  Elimination  of  the  Burke-Schumann  Restrictions 

The  restrictions  on  the  Burke-Schumann  analysis  prevent  its  application 
to  a  wide  range  of  problems.  In  this  section,  we  describe  how  the  computed 
results  are  affected  by  eliminating  some  of  the  restrictions  in  the  analysis. 
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In  these  computations,  the  ratio  of  the  inner  to  the  outer  radii  is  0.2  (a 
=  0.5,  b  =  2.5)  and  the  inlet  velocity  is  30  cm/s.  The  fuel  consists  of  3.41 
parts  of  #2  to  1  part  N2  by  volume  and  the  oxidizer  is  air. 

Introduction  of  correct  stoichiometry,  densities,  and  diffusion  coeffi¬ 
cients,  but  maintaining  the  conditions  of  uniform  inlet  velocity  and  isother¬ 
mal  reaction,  required  that  radial  gradients  be  included.  The  full  system 
of  equations  (1)  -  (4)  were  solved  for  this  case  and  the  results  are  shown 
in  Figure  4  after  0.4  s.  Figures  4a  and  4b  show  the  axial  and  radial  veloc¬ 
ity  contours  and  Figures  4c  and  4d  show  the  contours  of  fuel  and  oxidizer 
mole  fractions.  Small  radial  velocities  are  induced  in  the  reaction  zone  even 
though  heat  release  is  not  included.  The  maximum  axial  velocity  of  approx¬ 
imately  35  cm/s  occurs  about  0.25-0.50  cm  from  the  jet  centerline.  There 
is  a  region  of  lower  velocity  within  the  fuel-rich  zone  close  to  the  centerline 
which  initially  consisted  of  pure  fuel  mixture.  There  is  only  a  small  region 
very  close  to  the  inlet  which  is  still  pure  fuel  mixture  because  product  has 
diffused  across  the  jet.  This  diffusion  of  heavier  gases  increases  the  density 
in  the  fuel  zone.  Because  momentum  is  conserved,  the  velocity  decreases. 
The  flame  interface  lies  between  the  lowest  fuel  and  oxidizer  contours.  The 
flame  shape  is  slightly  distorted  due  to  the  fluctuations  in  the  radial  velocity. 

So  far,  heat  release  from  the  chemical  reaction  has  not  been  included, 
i.e.  Q  in  Equation  (3)  is  zero.  The  effects  of  including  heat  release  are 
shown  in  Figure  5  after  approximately  0.3  seconds.  The  volumetric  expansion 
associated  with  the  heat  release  accelerates  the  flow  resulting  in  maximum 
axial  velocities  of  120  cm/s.  The  hot  gases  are  accelerated  outward  in  the 
radial  direction  with  a  maximum  velocity  of  about  9  cm/s  (Fig.  5d).  The 
radial  velocity  contours  show  that  the  flame  is  not  steady  but  has  vortices 
which  form  at  the  fuel/oxidizer  interface.  This  flame  is  shorter  than  that 
predicted  when  heat  release  is  not  included  in  the  calculation.  The  increased 
radial  velocity  provides  am  additional  mechanism  for  mixing  so  the  reaction 
zone  is  wider  and  the  flame  is  shorter. 
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When  viscous  effects  are  included  in  the  calculation,  the  flow  field 
changes  again  and  these  results  axe  shown  in  Figure  6.  The  gradients  in 
the  axial  velocity  are  reduced  and  the  vortical  structures  are  damped.  The 
concentration  and  temperature  fields  axe  similar  to  those  without  viscosity. 
This  flame  is  slightly  longer  than  the  flame  without  viscosity  because  the 
radial  mixing  has  been  reduced. 

In  the  final  simulation,  gravitational  effects  were  included.  This  sim¬ 
ulation  includes  all  effects  discussed  in  the  equations  (1)  -  (4).  This  case 
was  started  from  the  flame  in  Figure  6  and  the  results  after  0.25  seconds  are 
shown  in  Figure  7.  The  radial  velocity  field  shows  large  structures  which 
form  in  the  high  temperature  region  near  the  jet.  These  structures  convect 
downstream  and  change  the  local  concentration  field.  The  H2  concentration 
field  is  not  affected  by  gravity  but  the  02  concentration  field  has  changed. 
Gravity  has  a  significant  effect  because  there  is  a  large  density  difference 
between  the  burnt  and  unburot  gases  in  this  flow. 

Figure  8  shows  a  time  sequence  of  02  mole  fraction  contours.  In  the  first 
frame,  a  bulge  appears  approximately  in  the  middle  of  the  computational 
domain  and  convects  upward  in  frames  2  and  3.  By  frame  4,  the  02  field  is 
very  distorted  as  the  structure  convects  upward.  In  the  final  two  frames,  the 
structure  continues  to  roll  up  as  it  convects  out  the  computational  domain. 
We  estimate  the  frequency  to  be  15-20  Hz.  These  structures  axe  similar  to 
those  observed  experimentally  in  unconfined  diffusion  flames  (20-21)  and  are 
often  attributed  to  the  effect  of  buoyancy. 

VI.  Conclusions 

The  Burke- Schumann  analysis  has  formed  many  of  our  ideas  about  lam¬ 
inar  diffusion  flames.  This  simplified  approach  describes  the  global  nature  of 
a  confined  laminar  flame  but  ignores  many  of  the  physical  phenomena  in  real 
flames.  In  this  paper,  we  described  a  new  computer  program  which  includes 
these  effects.  The  simulations  show  details  of  the  flames  which  cannot  be 
observed  from  the  analytical  solution. 
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Introducing  variable  density  and  diffusion  coefficients  for  a  Hi  —  Ni  fuel 
jet  with  coflowing  air  results  in  small  radial  velocities.  In  a  flame  where 
the  inlet  fuel  and  oxidizer  velocities  axe  equal,  heat  release  accelerates  the 
gases  and  produces  a  mixing  region  characterized  by  large-scale  instabilities 
which  are  damped  by  viscosity.  The  effects  of  heat  release  and  viscosity  are 
not  included  in  the  Burke- Schumann  analysis  but  they  appear  to  counter- act 
each  other. 

Gravity  produces  a  significant  change  in  the  flow  field  of  a  confined 
diffusion  flame.  The  flame  fluctuates  in  time  as  the  bouyancy-driven  struc¬ 
tures  convect  upward.  These  low  frequency  fluctuations  obviously  are  not 
represented  in  the  steady  state  analysis  of  Burke  and  Schumann,  but  these 
fluctuations  do  not  change  the  flame  location  significantly. 

Thus,  as  the  Burke-Schumann  restrictions  are  eliminated,  the  flame 
characteristics  change.  Realistic  stoichiometry,  diffusion  coefficients,  and 
densities  for  a  Hi  —  JV2  flame  with  coflowing  air  results  in  a  laminar  flame 
with  only  small  radial  velocities.  When  heat  release  is  included,  vortices  form 
in  the  reaction  zone  but  these  structures  are  damped  by  viscosity.  Finally, 
gravity  induces  large-scale  structures  to  form  in  the  region  outside  of  the 
reaction  zone. 
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Nomenclature 


Symbol 

Definition 

Speed  of  sound  (cm/s) 

c„ 

Specific  heat  (erg/g-K) 

Dik 

Binary  diffusion  coefficient  between  species  i  and  k  (cm2/s) 

E 

Total  energy  density  (erg/cm3) 

G 

Gravitational  acceleration  constant  (980.67  cm/s2) 

h 

Enthalpy  per  molecule  (erg/molecule) 

I 

Unit  tensor  (nondimensional) 

k 

Boltzmann  constant  (1.3805  x  10-16  erg/K) 

n 

Number  density  (cm-3) 

P 

Pressure  (dyne/cm2) 

Q 

Energy  release  rate  (erg/cm3  — s1 ) 

T 

Temperature  (K) 

V 

Velocity  (cm/s) 

w 

Production  rate  (cm-3s-1) 

X 

Spatial  coordinate  (cm) 

X 

Mole  fraction 

y 

Spatial  coordinate  (cm) 

Y 

Mass  fraction 

Greek 

e 

Specific  internal  energy  (erg/cm3) 

7 

Ratio  of  specific  heats,  cp/cv 

K 

Thermal  conductivity  coefficient  (erg/s-K-cm) 

P 

Coefficient  of  shear  viscosity  (poise,  g/cm-s) 

P 

Mass  density  (g/cm3  ) 

T 

Viscous  stress  tensor  (dynes /cm2) 

Superscripts 

T 

Transpose  operation  on  a  matrix 

Subscripts 

m 

Mixture  of  species 

i,j ,  fc,  or  / 

Individual  species 
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Fig.  1  —  Geometry  used  for  Burke-Schumann  calculations  showing  flame  location  for  typical  underventilated  or  overven¬ 
tilated  flame 
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—  Contours  of  fuel  concentration  normalized  by  inlet  fuel  concentration  for  overventilated  Burke-Schumann  at  times 
(b)  0. 1  sec  (c)  0.4  (e)  1.0  seconds 
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Fig.  2b  —  Contours  of  oxidizer  concentration  normalized  by  inlet  oxidizer  concentration  for  overventilated  Burke-Schumann 
at  times  (a)  0.0  (b)  0. 1  sec  (c)  0.4  sec  (d)  0.7  (e)  1 .0  seconds 
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Pig.  3  —  Comparison  of  analytical  and  computed  solutions  for  flame  location  for  three  different  Burke-Schumann  flames, 
(a)  Pure  fuel  reacting  with  1  part  oxidizer  +  1  part  inert,  (b)  One  part  fuel  +  3  parts  inert  reacting  with  pure  oxidizer,  (c) 
Pure  fuel  reacting  with  1  part  oxidizer  +  4  parts  inert. 
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g  7  _  Contours  of  (a)  radial  velocity  (b)  axial  velocity  (c)  mole  fraction  H2  (d)  mole  fraction  02  (e)  temperature  for 
2  -  N2  diffusion  flame  with  heat  release,  viscosity,  and  gravity.  Dimensions  are  in  cm,  velocities  are  in  cm/s,  tempera¬ 
ture  is  in  K. 
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Fig.  8  —  Time  sequence  for  contours  of  02  mole  fraction  showing  the  formation  of  large  structure.  The  time  interval 
between  each  frame  is  10  ms. 


